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Abstract 
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I. INTRODUCTION 



In hadron-hadron collisions, the transverse momentum of Drell-Yan pairs or weak gauge 
bosons (W ± and Z) is generated by emission of gluons and quarks, as predicted by Quan- 
tum Chromodynamics (QCD). Therefore, in order to test the QCD theory or electroweak 
properties of the vector bosons, it is necessary to include the effects of multiple gluon emis- 
sion. One theoretical framework designed to account for these effects is the resummation 
formalism developed by Collins, Soper, and Sterman (CSS) [1], which has been applied to 
study production of single [2-5, and references therein] and double [6, 7] electroweak gauge 
bosons, as well as Higgs bosons [8], at hadron colliders. Just as the nonperturbative func- 
tions, i.e., parton distribution functions (PDF's), are needed in order to predict inclusive 
rates, an additional nonperturbative function, W^, p , is required in the CSS formalism to 
describe the transverse momentum of, say, weak bosons. Many studies have been performed 
in the literature to determine W^ p using the available low energy Drell-Yan data [5, 9-13]. 
In particular, in Ref. [11] three of us have examined various functional forms of W^ p to 
test the universality of the CSS formalism in describing the Drell-Yan and weak boson data. 
The result of that study was summarized in Table II of Ref. [11]. In addition, that paper 
made several important observations, as to be discussed below. First, neither the DWS form 
(cf. Eq. (10)) nor the LY form (cf. Eq. (11)) of the nonperturbative function W^ p could 
simultaneously describe the Drell-Yan data in a straightforward global fit of the experiments 
R209 [14], E605 [15], and E288 [16], as well as the Tevatron Run Z boson data from the 
CDF collaboration [17]. Hence, it was decided in Ref. [11] to first fit only the first two 
mass bins (7 < Q < 8 GeV and 8 < Q < 9 GeV) of the E605 data and all of the R209 
and the CDF-Z boson data in the initial fits A 2 and A 3 . In total, 31 data points were 
considered. We allowed the normalization of the R209 and E605 data to float within their 
overall systematic normalization errors, while fixing the normalization of the CDF-Z Run 
data to unity. Second, after the initial fits, we calculated the remaining three high-mass 
bins of the E605 data not used in Fits and found a reasonable agreement with the 
experimental data. In order to compare with the E288 data, we created Fits A^ 2 3 , in which 
we fixed the functions W%, p to those obtained from Fit A 3 and performed a fit for NORM 
(the fitted normalization factor applied to the prediction curves for a given data set) for the 
E288 data alone. We found that the quality of the fit for the E288 data is very similar to 
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that for the E605 data, and that the normalizations were then acceptably within the range 
quoted by the experiment. Hence, we concluded that the fitted functions W^ p reasonably 
describe the wide-ranging, complete set of data, in the sense discussed above. Third, most 
importantly, we found that the complete set of data available in that fit was not yet precise 
enough to clearly separate the g 2 and g\g^ parameters within a pure Gaussian form of W^ p 
with x dependence similar to that of the LY form. This Gaussian form is given explicitly in 
Eq. (12). 

In this paper, we show that, after including the transverse momentum distributions of the 
Z bosons measured by the D0 [18] and CDF collaborations [19] in Run-1 at the Fermilab 
Tevatron, we are able for the first time to perform a truly global fit of the nonperturbative 
function W^ p to the complete set of data on vector boson production. In this fit, the data 
from the experiments R209, E605, and E288, as well as the Tevatron Run-1 Z data are 
treated on the same footing. We emphasize that in this new fit, the E288 data are also 
included in the global fit, in contrast to the study done in Ref. [11]. Furthermore, we show 
that the Gaussian form of W^ p given in Eq. (12) clearly fits the data the best, as compared 
to either the updated DWS form (10) or updated LY form (11). These nice features are 
driven by the inclusion of the Run 1 Z data, for these data determine the g 2 coefficient with 
good accuracy by separating the contributions from g 2 and gig% . 

This paper is organized as follows. In Section II, we briefly review the CSS formalism 
with an emphasis on its nonperturbative sector. In Section III, we describe in detail the 
results of our fits. In Section IV, we discuss various aspects of our study and comment on the 
validity of our approach to the treatment of the nonperturbative region. Section V contains 
conclusions. 



II. COLLINS-SOPER-STERMAN RESUMMATION FORMALISM 

As an example, we consider production of a vector boson V in the collision of two hadrons 
hi and h 2 . In the CSS resummation formalism, the cross section for this process is written 
in the form 



da(h!h 2 -> VX) 1 



dQ 2 dQ 2 T dy (2vr) 



; 5(Q 2 -M*) j d 2 be^ Ti; W fk (b,Q,x u x 2 ) + Y(Q T ,Q,x u x 2 ), (1) 



2 



where Q, Qt, and y are the invariant mass, transverse momentum, and the rapidity of the 
vector boson V. The Born-level parton momentum fractions are defined as x\ = e v Q/\J~S 
and x 2 = e~ y Q/\/S, with \/S being the center-of-mass (CM) energy of the hadrons h\ and 
h 2 . We will refer to the integral over the impact parameter b in Eq. (1) as the "W-term". Y is 
the regular piece, which can be obtained by subtracting the singular terms from the exact 
fixed-order result. The quantity Wj% satisfies a set of renormalization- and gauge- group 
equations [20] with the solution of the form 

W fk (b, Q, Xl , x 2 ) = e~ s ^ b ' c ^W fk (b, |l Xl , x 2 \ , (2) 



where C\ and C 2 are constants of order unity, and the Sudakov exponent is defined as 



S{Q,b,C u C 2 ) 



c ^ 2 dff 



A {a s {ji), C x ) In ( ) + B (a s {]l), C Xl C 2 ) 



(3) 



The dependence of (b, ■^,x 1 ,x 2 S j on x\ and x 2 factorizes as 

W fk (b, x x , x 2 j = f Ui, b, W- kh2 (x 2 , b, + (j <-> k). (4) 

In the perturbative region, i.e., at b 2 <ti 1/Aq CD , the function 3? jh can be expressed as 
a convolution of the parton distribution functions f a /h with calculable Wilson coefficient 
functions Cj a : 

w * ( X > b ^) = T^lj c - (f > h %^ = t) fa/h ( e ' ^ = t) • (5) 

The sum over the index a is over all types of incoming partons. The sum over the index 
j(k) is over all quarks (antiquarks) . The coefficient o"o includes constant factors and quark 
couplings from the leading-order cross section, which can be found, e.g., in Ref. [1]. The 
factorization scale \i on the r.h.s. of Eq. (5) is fixed to be C%/b. 
A few comments should be made about this formalism: 

• If both Q and 1/6 are much larger than the typical internal hadronic scale Aqcd, 
the A, B and C functions can be calculated order- by-order in a s . In our fit, we shall 
include the A and B functions up to 0(a 2 s ), and C functions up to 0(a s ). 

• A special choice can be made for the renormalization constants C{ to remove some of 
the logarithms in Wjk(b, Q, x\, x 2 ). This canonical choice is C\ = C3 = 2e _7E = 60 



and C 2 = C1/&0 = 1, where 7^ is the Euler's constant. We shall use this choice in our 
calculations. 

• In Eq. (1), the variable b is integrated from to 00. When b > 1 GeV -1 , the per- 
turbative calculation for W^Q), Q, Xi,x 2 ) is no longer reliable, and complicated long- 
distance physics comes in. Furthermore, even in the perturbative region (b < 1 GeV" 1 ) 
Wjk(b } Q, Xi, x 2 ) may contain some small nonperturbative terms, which arise, e.g., from 
power corrections to the CSS evolution equations. It is important to emphasize that 
significance of nonperturbative contributions of both types is drastically reduced when 
Q is of order of the W, Z boson masses or higher [1]. For those Q, the most part of the 
Qt distribution can be predicted based purely on the perturbative calculation, with 
the exception of the region of Qt below a few GeV, where sizeable sensitivity to the 
nonperturbative input remains [2, 9]. 

According to the common assumption, nonperturbative contributions to Wj%(b, Q, x±, x 2 ) 
can be approximated by some phenomenological model with measurable and universal 1 pa- 
rameters. Collins, Soper, and Sterman [1] suggested the introduction of the nonperturbative 
terms in the form of an additional factor W^ p (b, Q,Xi,x 2 ), usually referred to as the "non- 
perturbative Sudakov function". More precisely, the form factor Wji(b, Q, x 1; x 2 ) in Eq. (1) 
is expressed in terms of its perturbative part Wj^ rt and nonperturbative function W^ p as 2 

= ^;r( & *)^rw , (6) 

with 

K = b (7) 

+ (b/b max y 

In numerical calculations, b max is typically set to be of order of 1 GeV -1 . The variable 
b* never exceeds b max , so that W^ rt (b*) can be reliably calculated in perturbation theory 
for all values of b. Based upon the renormalization group analysis, Ref. [1] found that the 
nonperturbative function can be generally written as 



Wjl (b, Q, Q , xi, x 2 ) = exp 



In - F j/hl (x lt b) - F~ k/h2 (x 2 , b) 



(8) 



1 Here, we mean "universal" in the context of Drell-Yan-like processes, in which the initial state of the Born 
level process involves only quarks and antiquarks, and the observed final state does not participate in 
strong interactions. 
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Here and after, we suppress the arguments Q, X\ and x%, and denote W^(b,Q,x±,X2) as W^(b), etc. 
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where F\, Fj/h and Fu/u must be extracted from the data, with the constraint that 



^(6 = 0) = !. 



(9) 



Furthermore, F\ depends only on b. Fj/ hl and Fy h2 in general depend on x\ or x 2 , and their 
values can depend on the flavor of the initial-state partons (j and k in this case). Later, 
it was shown in Ref. [21] that the Fi(6)ln dependence is also suggested by infrared 

renormalon contributions to the form factor Wjj.(b, Q, Xi, x 2 ). The CSS resummation for- 
malism suggests that the nonperturbative function is universal. Its role is analogous to that 
of the parton distribution function in any fixed-order perturbative calculation. In particular, 
its origin is due to the long-distance effects that are incalculable at the present time, and 
its value must be determined from data. 

As discussed in Ref. [11], we will consider three different functional forms for W^ p . 
They are 



the 2-parameter pure Gaussian form, called the Davies- Webber-Stirling (DWS) form 
[9], 



exp 



-gi - 92 In 



Q_ 

2Q 



(10) 



the Ladinsky-Yuan (LY) form [10], 



exp 



-gi - g 2 In 



Q_ 

2Qc 



6 2 -[^ 3 ln(100x i; r 2 )]6 



which has a logarithmic x-dependent term linear in b; 



(11) 



and the 3-parameter pure Gaussian form, called the Brock-Landry-Nadolsky-Yuan 
(BLNY) form, 

Q 



exp 



-gi - g 2 In 



2Qc 



5-15-3 In (lOOxi^) 



b 2 . 



(12) 



We will refer to the updated DWS and LY parameterizations obtained in the current global 
fit as "DWS-G" and "LY-G" parameterizations, respectively, to distinguish them from the 
original DWS and LY parameterizations [9, 10] obtained in (non- global) fits to a part of the 
current data. 
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Table I: Vector boson production data used in this analysis. Here, SN exp is the published normal- 
ization uncertainty for each experiment. 



Experiment 


Reference 


Reaction 


VS (GeV) 


SN exp 


R209 


[14] 


P + p + x 


62 


10% 


E605 


[151 

L J 


p + Cu^ u + u~ + X 


38.8 


15% 


E288 


[16] 


p + Cu^ + x 


27.4 


25% 


CDF-Z 
(Run-0) 


[171 

L J 


p + p^ Z + X 


1800 




D0 -Z 
(Run-1) 


[18] 


p + p^ z + x 


1800 


4.3% 


CDF-Z 
(Run-1) 


[19] 


p + p^ z + x 


1800 


3.9% 



III. RESULTS OF THE GLOBAL FITS 

In order to examine the impact of including the Z boson data from the Run-1 of the 
Tevatron on the global fits and compare the new results to those given in Ref. [11], our 
theory calculations will consistently use CTEQ3M PDF's [22]. 3 For the same reason, we 
take Qo = 1.6 GeV and b max = 0.5 GeV -1 in all fits. 

As discussed in the previous sections, our primary goal is to determine the nonper- 
turbative function of the CSS resummation formalism. Hence, we need to include those 
experimental data, for which the nonperturbative part dominates the transverse momen- 
tum distributions. This requirement suggests using the low-energy fixed-target or collider 
Drell-Yan data in the region where the transverse momentum Qt of the lepton pair is much 
smaller than its invariant mass Q. Because the CSS formalism better describes production 
of Drell-Yan pairs in the central rapidity region (as defined in the center-of-mass frame of 
the initial-state hadrons), we shall concentrate on the data with those properties. Based 
upon the above criteria, we chose to consider data from the experiments listed in Table I 
and in kinematical ranges shown in Table II. We have also examined the E772 data [23] 

3 In principle, the non-perturbative function depends on the choice of the PDF's. However, we will argue 
later that this dependence can be currently neglected within the accuracy of the existing data. 



Table II: The data sets used for the fit. Pt and Q denote the published transverse momentum and 
mass of the Drell-Yan pair or the Z boson, respectively. 



Experiment 


Pt range 


Q range 




(GeV) 


(GeV) 


R209 


0.0 - 1.8 


5.0 - 11.0 


E605 


0.0 - 1.4 


7.0 - 9.0 & 10.5 - 18.0 


E288 


0.0 - 1.4 


5.0 - 9.0 


CDF-Z 


0.0 - 22.8 


91.19 


(Run-0) 






D0 -Z 


0.0 - 22.0 


91.19 


(Run-1) 






CDF-Z 


0.0 - 22.0 


91.19 


(Run-1) 







from the process p + H 2 — ► +X at v^S" = 56.6 GeV and found them incompatible with 
the rest of the data sets. Hence the E772 data were not included in the presented fits. 4 

The theoretical cross sections were calculated with the help of the resummation package 
Legacy, which was also used in previous fitting [10, 11] and analytical studies [3, 4, 6, 7, 25- 
27], as well as for generating input cross section grids for ResBos Monte-Carlo integration 
program [4]. This package is a high-performance tool for calculation of the resummed cross 
sections, with the computational speed increased by up to a factor 800 after the reorganiza- 
tion and translation of the source code into C/C++ programming language in 1999-2001. 
During the preparation of this paper, we confirmed the stability of the numerical calcula- 
tion of the resummed cross sections (1) by comparing the output of several Fourier-Bessel 
transform routines based on different algorithms (adaptive integration, Fast Fourier-Bessel 
transform [28], and Wolfram Research Mathematica 4.1 NIntegrate function). Specifi- 
cally, the outputs of three routines are in a very good agreement at all values of Qt- For 
instance, the Z boson cross sections presented in this paper and Ref. [11] are calculated with 

4 For the best-fit values of gi given below, the theory prediction for the E772 experiment is typically smaller 
than the data by a factor of 2. Similarly, the E772 data are not well fit in the CTEQ global analysis of 
parton distribution functions [24]. 
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the relative numerical error less than 0.5% at Qt < 50 GeV and less than 1 — 2% at Qt > 50 
GeV. Note that the relative error of about 1% is comparable with the size of higher-order 
(NNLO) corrections, as well as numerical uncertainties in the existing two-loop PDF sets. 
More details on the tests of accuracy of the resummation package will be presented elsewhere 
[29]. 5 

Using the above sets of the experimental data, we fit the values of the nonperturbative 
parameters g\, gi and in the DWS-G form (10), LY-G form (11), and BLNY form (12) 
of the nonperturbative function W^ p (b, Q, Qo, x±, x 2 ). Since we allow the normalizations 
for the data to float within the overall systematic normalization errors published by the 
experiments, the best-fit values of g±, g 2 and g% are correlated with the best-fit values of 
the data normalization factors Nf it (individually applied to each data set). Note that the 
normalization of the CDF-Z Run data was fixed to unity due to their poor statistics as 
compared to the Run-1 data. 

Table III summarizes our results. To illustrate the quality of each fit, Figs. 1-5 compare 
theory calculations for the DWS-G, LY-G, and BLNY parameterizations to each data set. 
We emphasize again that the new LY-G parameterization presented in Table III was obtained 
by applying the conventional global fitting procedure to the enlarged data set listed in 
Tables I and II. In contrast, the original LY fit in Ref. [10] was obtained by first fitting the 
g 2 parameter using the CDF-Z Run and R209 data, and then fitting g\ and g 3 parameters 
after including the other available Drell-Yan data (which is a small subset of the data in the 
current fit). 

It is evident that the Gaussian BLNY parameterization fits the whole data sample no- 
ticeably better than the other two parameterizations, both in terms of \ 2 P er degree of 
freedom (dof) given in Table III and in terms of the pictorial comparison in Figs. 1-5. When 
compared to the Run-1 Z data, the DWS-G and LY-G fits both fail to match the height and 
position of the peak in the transverse momentum distribution (Figs. 4 and 5). Similarly, 
according to Fig. 3, the BLNY parameterization leads to the best agreement with the E288 
data. The three fits are indistinguishable when compared to the E605 data, cf. Fig. 2. Fig. 1 
shows a clear difference between the BLNY fit and the other two fits for the lowest mass bin 

5 An interface to the simplified version of Legacy and online plotter of resummed transverse momentum 
distributions are available on the Internet at http://hep.pa.msu.edu/wwwlegacy/ . 



Table III: The results of the fits. Here, Nf it is the fitted normalization for each experiment. 



Parameter 


DWS-G fit 


LY-G fit 


BLNY fit 


91 


0.016 


0.02 


0.21 


92 


0.54 


0.55 


0.68 


93 


0.00 


-1.50 


-0.60 


CDF Z Run-0 


1.00 


1.00 


1.00 


Ntit 


(fixed) 


(fixed) 


(fixed) 


R209 


1.02 


1.01 


0.86 


Ntu 
i "fit 








E605 


1.15 


1.07 


1.00 










E288 


1.23 


1.28 


1.19 


Nut 








D0 Z Run-1 


1.01 


1.01 


1.00 


Nf U 








CDF Z Run-1 


0.89 


0.90 


0.89 


N fa 








x 2 


416 


407 


176 


X 2 /dof 


3.47 


3.42 


1.48 



a Thus, by definition, NORM in Ref. [11] is equal to 1/N flt . 

of the R209 data (the upper data in the Figure). However, the x 2 contribution from this 
mass bin is about the same for all three fits. 

The error on the fitted nonperturbative parameters g±, gi and g% can be calculated by- 
examining the x 2 distribution of the fit. For the BLNY form 

Q 



W jl P ( b i Qi Qo, xi,x 2 ) = exp 



2Q 



(13) 



-^(l + ^M 1003 ^)) ~92 In 
with Q = 1.6 GeV and b max = 0.5 GeV~ , we found that 

9l = 0.21±Hi GeV 2 , g 2 = 0.68t°;°2 GeV 2 , £ 3 = -0.6^. 
The errors in Eq. (14) were computed as follows. First, x 2 values were calculated around their 
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(14) 



R209 Data 

140 | , . , , , , , . , , . , ^- 




' 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 • 1 1 

0.5 1 1.5 2 

P y (GeV) 

Figure 1: Comparison to the R209 data for the process p+p — ► fi + fi~ + X process at y/S = 62 GeV. 
The data are the published experimental values. The curves are the results of the fits and are 
multiplied by the best-fit values of 1/Nfn given in Table III. 

minimum Xmin m order to obtain, in essence, a three dimensional function x 2 (.9i, 9i, <?3)- 6 
Next, we plotted an ellipsoid surface determined by the condition X 2 (9i, 9i, #3) = Xmin + 1 
in the three-dimensional space of parameters gi, cf. Fig. 6. The extremes in each coordinate 
for this surface were taken as the errors for the respective parameters. Finally, we note that 
using xLin + 1 as the confidence limit for determining the values of gi in the presence of sub- 
stantial systematic errors is generally idealistic, for the experiments often make judgments 
on systematic uncertainties that are not of a Gaussian nature. Further discussion of this 
issue can be found in the next Section. 

6 In this calculation, we scanned the values of gi and 172 between and 1, and 53 between —2 and 3. 
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E605 Data 



o 
II 



0.1 



T3 | "5r 



0.01 




Data 

Normalized LY-G Fit 
Normalized DWS-G Fit 
Normalized BLNY Fit 




1 .4 



Figure 2: Comparison to the E605 data for the process p + Cu — > /x + pC + X at y/S = 38.8 GeV. 
The data are the published experimental values. The curves are the results of the fits multiplied 
by the best-fit values of 1/Nfn given in Table III. 

IV. DISCUSSION 

To fit the complete set of the experimental data (more than 100 data points), this analysis 
introduced 3 free parameters (gi, g 2 and g 3 ) in the parameterization of W^ p , together with 
the chosen values of the parameters Q and b max . In this paper and Ref. [11], we have chosen 
Qo = 1.6 GeV, which coincides with the lowest energy scale used in the CTEQ3M PDF set. 
This choice is nothing more than a matter of convenience, since it enforces positivity of 
the logarithm \n(Q/2Q ) in the Q range of the fitted data (Q > 5 GeV). But Qo does 
not have to take that special value, since its variations can be compensated for in the full 
nonperturbative function W^ p by adjusting parameters gi and g 3 (at g 2 fixed). For instance, 
we could have chosen Q to be equal to l/b max and rewrite Eq. (12) as 



Wjl P (b,Q,Q ,x 1 ,x 2 ) =exp -^(1 + g' 3 In (100xix 2 )) - g 2 ln(Q6 n 



(15) 
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E288 Data 



Data 

Normalized LY-G Fit 

— — Normalized DWS-G Fit 
Normalized BLNY Fit 



m 
II 



> 
o 



J= he, 

"a I -a* 



0.1 




0.2 



0.4 



0.6 



0.8 



1 .2 



1 .4 



P y (GeV) 



Figure 3: Comparison to the E288 data for the process p + Cu — > fi + fi~+X at y/S = 27 .4 GeV. The 
data are the published experimental values. The curves are the results of the fits and are multiplied 
by the best-fit values of 1/Nfn given in Table III. 



where g 2 = g 2 , g[ = g\ - g 2 \n(2Q b max ), and g' 3 = g^gxj g\. This new form of W^ p would be 
completely identical to the original form in Eq. (12). Hence, the total number of parameters 
needed to describe Wji{b) in the used prescription is four, i.e., gi and b max . 

While the parameter Q plays no dynamical role, the meaning of the parameter b max is 
quite different. Roughly speaking, its purpose is to separate nonperturbative effects from 
perturbative contributions through the introduction of the variable b* defined in Eq. (7). 
According to its definition, the variable 6* is practically equal to b when b 2 -C b 2 max . For 
b — ► oo, it asymptotically approaches b max . Hence, in the factorized CSS representation (6) 
the perturbative part W^ rt (b*) approaches its exact value (evaluated at b) when b — > 0, and 



it is frozen at 6* 



when b — > oo. While b max should lie in the perturbative region to 



make the computation of WJ£ T (6*) feasible, it is also desirable to make it as large as possible, 
to reduce deviations of W P j Tt from its exact behavior at smaller b. Note, however, that the 
changes in the behavior of W p ? rt due to the 6* prescription, such as the loss of suppression of 
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1 5 




20 



Figure 4: Comparison to the D0 -Z Run-1 data. The data are the published experimental values. 
The curves are the results of the fits and are multiplied by the best-fit values of 1/Nfa given in 
Table III. 

Wj^ rt by the perturbative Sudakov factor at b > b max , can be compensated for by increasing 
the magnitude of the nonperturbative function W^ p . Hence, in the 6* prescription the 
nonperturbative function W^ p generally serves a dual purpose of the parameterization for 
truly nonperturbative effects and compensation factor for modifications in Wj^ rt due to the 
6* variable. Consequently, the best-fit parameterization of W^ p depends on the choice of 

bmax ■ 

Based on the fact that the 6* prescription with the BLNY form of the non-perturbative 
function provides an excellent fit to the whole set of Drell-Yan-like data, 7 we conclude 
that the 6* prescription remains an adequate method for studies of resummed Qt distribu- 



7 For example, as shown in Figs. 4 and 5, both the peak position and the shape of the transverse momentum 
distribution of the Z boson measured by the D0 and CDF Collaborations at the Run-1 of the Tevatron 
are very well described by the theory calculations. This feature is extremely important for the precision 
measurement of the W boson mass at the Tevatron. 
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CDF Z Run 1 

800 I 1 1 1 1 1 1 1 1 1 1 1 1 1 r 




' — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 
5 10 15 20 

P t (GeV) 

Figure 5: Comparison to the CDF-Z Run-1 data. The data are the published experimental values. 
The curves are the results of the fits and are multiplied by the best-fit value of 1/Nfu given in 
Table III. 

tions.This, however, does not mean that future tests and improvements of this method will 
not be necessary. For instance, one might consider a larger value of b max for the fit. Up to 
now, all the fits performed in the literature have taken b max to be 0.5 GeV -1 . However, b max 
could be chosen to be as large as b /jip, where b = 2e -7E ~ 1.122..., and /Jp is the initial 
energy scale for the PDF set in use. Several of the most recent PDF sets, such as CTEQ6M 
[32], MRST'2001 [33] and MRST'2002 [34], use fi° F as low as 1 GeV. Correspondingly, for 
the new PDF sets b max can be as large as 1.122 GeV -1 . It will be interesting to test the 
6* prescription with such increased value of b max in a future study. Moreover, an approach 
alternative to the fe* prescription was recently proposed in Refs. [12, 13]. Fig. 7 of Ref. [13] 
shows the comparison of this new theory calculation to the Run-1 CDF data, which should 
be compared to Fig. 5 of this paper. Furthermore, in Refs. [30, 31], another method for per- 
forming the extrapolation to the large-6 region was proposed. To see the difference between 
the theory predictions for the Tevatron Z data, one can compare Figs. 5 and 6 of Ref. [31] 



14 



. 7 



g2 
0.68 



0.66 



-0.55 




Figure 6: Uncertainty contour and two dimensional projections for the 3-parameter Gaussian BLNY 
form fit, cf. Eq. (12). 
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Figure 7: Transverse momentum distributions of Z bosons at the Tevatron Run-1. The theory 
curves, calculated using the CTEQ6M PDF's [32] and the BLNY parameterization (solid line) or 
the original LY parameterization, shown as dashed line, are compared to the D0 data and CDF 
data. The data curves have been multiplied by an overall normalization factor 1.0 in D0 data and 
0.9 in CDF data. 
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to Figs. 5 and 4 of this paper. In a forthcoming publication [29], we shall present a detailed 
comparison of various methods for description of large-6 physics to the Tevatron Z data. 

We would like to conclude this section with a remark about the dependence of our results 
on the choice of parton distribution functions. As noted above, the fitted nonperturbative 
function W^ p (b) in the CSS resummation formalism is correlated to the PDF's used in the 
theory calculation. In the current fit, we have chosen to use CTEQ3M PDF's to facilitate 
the comparison with the previous results in Ref. [11]. The usage of a more modern PDF set 
with the BLNY parameterization of W^ p would result in a difference of a few percent in 
Qt distributions. Currently, such differences are of order of normalization errors quoted by 
the experiments. Hence, good agreement between the theory and data can be obtained by 
small normalization shifts in the data. We illustrate this point in Fig. 7, which compares 
the Tevatron Run-1 Z data to the CSS resummation calculation performed using CTEQ6M 
PDF's [32] and BLNY parameterization (12) with the best-fit nonperturbative parameters 
(14). By adjusting the normalization of the CDF data by the best-best value Nf it = 0.89 
in Fig. 7, the theory is brought in a good agreement with both sets of the Run-1 Z data. 
For comparison, we also show in Fig. 7 the prediction from using CTEQ6M PDF's and the 
original LY parametrization [10]. 

As the quality of the data improves in the future, the correlation between the nonpertu- 
bative function W^ p {b) and PDF's will become more important. Hence, it will be certainly 
desirable to repeat the work done in this paper using the most recent set of PDF's and 
potentially perform the joint global analysis of the PDF's and CSS nonperturbative func- 
tion. Furthermore, as the statistical errors decrease, correct treatment of the systematic 
errors (which are hardly of a Gaussian nature in many experiments) becomes ever more 
necessary. Therefore, a method more elaborate than the simple criterion Xmin + 1 used 
in the current analysis will be needed to determine the true confidence limit for the non- 
perturbative function. Recently, new efficient methods were proposed to perform the error 
analysis on the nonperturbative parameters in the PDF's in the presence of systematic errors 
[32, 34, 35]. In the future, the same methods can also be applied to determine the errors for 
the nonperturbative parameters gi, gi y and g% in the BLNY nonperturbative function. 
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V. CONCLUSIONS 



We have shown for the first time that the complete set of low energy Drell-Yan data 
(R209, E605 and E288) and the Tevatron Run-1 Z-boson data can be simultaneously de- 
scribed by the CSS resummation formalism. This is the first truly global fit, which treats 
all the low energy Drell-Yan data and the Tevatron Z data on the same footing. This fact 
strongly supports the universality of the CSS nonperturbative function W^ p . Just as the 
universality of the PDF's allows us to predict the inclusive rate for a scattering process in 
hadron collisions, the universality of W^ p allows us to predict the transverse momentum 
distributions in production of single vector bosons and other similar processes at hadron 
colliders. For example, if true, the resummation formalism can be applied not only to the 
Drell-Yan pair production and W (or Z) boson production, but also to associated produc- 
tion of Higgs and W (or Higgs and Z) bosons [26], diphoton production [6], ZZ (or WW) 
pair production [7], and s-channel neutral or charged Higgs boson production (induced from 
quark fusion with large Yukawa coupling to Higgs boson) [27]. When these processes are 
measured to a good accuracy, they can be used to further test the CSS resummation for- 
malism. Just as the PDF's must be constantly refined in order to fit new experimental 
data, the nonperturbative function W^, F may also require modification in the future. At 
the current stage, the existing data shows remarkable preference for a Gaussian form of the 
nonperturbative function with Q dependence predicted by the CSS formalism. 
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